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Abstract. The quasistatic rate-independent damage combined with linearized plasticity with hard¬ 
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1 INTRODUCTION 

A combination of plasticity and damage, also called ductile damage , open colorful scenarios 
with important applications in civil or mechanical engineering and with interesting mathe¬ 
matical problems, in particular in comparison with mere plasticity or mere damage. Often, 
both plastihcation and damage processes are much faster than the rate of applied load and, 
in a basic scenario, any internal time scale is neglected and the mentioned inelastic processes 
are considered as rate independent. The goal of this article is to devise a model together 
with its efficient computational approximation that would lead to a numerical stable and 
convergent scheme and, at least in particular situations, calculate a physically relevant so¬ 
lutions of a stress-driven type verifiable aposteriori by checking a suitable version of the 
maximum-dissipation principle. 

We use the very standard linearized, associative, plasticity at small strain as presented 
e.g. in pT]. Simultaneously, we use also a rather standard scalar (i.e. isotropic) damage as 
presented e.g. in uni. We have primarily in mind a conventional engineering model with 
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unidirectional evolution of damage; actually, the healing will be here allowed rather from 
analytical reasons and can be expected ineffective in usual applications, cf. Remark 12.11 
below. All rate-dependent phenomena (as inertia or heat conduction and thermo-coupling) 
are neglected, this means the problem is considered as quasistatic and fully rate-independent. 
To avoid serious mathematical and computation difficulties, we have in mind an incomplete 
damage. 

The mentioned modelling simplification leading to quasistatic rate-independent system 
however, which reflects a certain well-motivated asymptotics, brings however quite serious 
questions and difficulties because the class of reasonably general solutions is very wide if the 
governing energy is not convex (as necessarily here) and involves solutions of very different 
nature, some of them physically not relevant; cf. [20, [21]. In particular, to avoid unwanted 
effects of unphysically easy damage under subcritical stress, one cannot require (otherwise 
attractive idea of) energy conservation and thus cannot consider so-called energetic solutions 
in the sense [22], although this concept is occasionally used for damage with plasticity in 
purely mathematically-focused literature mm- This is related with the discussion whether 
rather energy or rather stress is responsible for governing evolution of rate-independent 
systems [!§] . 

In contrast to the mentioned energetic solutions (which allows for simpler analysis with¬ 
out considering gradient plasticity but lead to recursive global-minimization problems which 
are difficult to realize and may slide to unphysically scenarios of unrealistic early damage), 
we will focus here on solutions that are rather stress driven and that can be efficiently ob¬ 
tained numerically. We will rely on a certain careful usage of a suitable integral-version of the 
maximum-dissipation principle, as devised in [31] and used, rather heuristically, in engineer¬ 
ing models of damage with plasticity and hardening, cf. [8]. This brings specific difficulties 
with convergence (which requires usage of gradient plasticity) and specific a-posteriori veri¬ 
fication of a suitable approximate version of the mentioned maximum-dissipation principle, 
as suggested in [21] Rem. 4.6] for damage itself, modified here for the combination of dam¬ 
age and plasticity analogously like in [35] for a surface variant of the elasto-plasto-damage 
model. If the maximum-dissipation principle holds (at least with a good accuracy) we can 
claim that the numerically obtained solution is physically relevant as stress-driven (with a 
good accuracy). 

A physically more justified and better motivated approach would be to involve a small 
viscosity to the damage variable or to the elastic and the plastic strains, and then to pass 
these viscosities to zero. The limits obtained by this way are called vanishing-viscosity 
solutions to the original rate-independent system and their analysis and computer imple¬ 
mentation is very difficult; for models without plasticity cf. mi for viscosity in damage or 
[34] for viscosity in elastic strain. 

In principle, there are two basic scenarios how the material might respond to an increasing 
loading: either first plasticize and then go into damage due to hardening effects , or first go 
into damage and then plasticize; of course, various compromising scenarios are possible too. 
The latter scenario needs a damage influenced yield stress and allow for no hardening (and in 
particular perfect plasticity), cf. [35] . Let us only remark that a damage-dependence of the 
yield stress in the fully rate-independent setting would make the dissipation state-dependent, 
which brings serious difficulties as seen e.g. in [2TJ Sect. 3.2] and, for the particular elasto- 
plasto-damage model, in p, [3] id]. In this paper, we will however concern exclusively on the 
former scenario, i.e. in particular, the damage does not influence the yield stress. Moreover, 
we will consider only kinematic hardening, although all the considerations could easily be 
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augmented by isotropic hardening, too. (Another essential difference from [36] is that, as 
already explained, the energy is intentionally not conserved in here.) 

The plan of the paper is as follows: In Section [2] we devise the model in its classical for¬ 
mulation and then, in Section [3j its suitable weak formulation with discussing stress-driven 
solutions and the role of the maximum-dissipation principle. In Section Qj we propose a 
constructive time discretisation method and prove its numerical stability (i.e. a-priori esti¬ 
mates) and convergence towards weak solutions. After a further fini te-element, discretisation 
outlined in Section [5j this allows for efficient computer implementation of the model, which 
is demonstrated on illustrative 2-dimensional examples in Section [6j 


2 THE MODEL AND ITS WEAK FORMULATION 


Hereafter, we suppose that the damageable clasto-plastic body occupies a bounded Lipschitz 
domain 12 C M d , d = 2 or 3. We denote by ft the outward unit normal to <912. We further 
suppose that the boundary of 12 splits as 

<912 ■= r = r D u r N , 


with T d and T N open subsets in the relative topology of <912, disjoint one from each other and, 
up to {d— l)-dimensional zero measure, covering <912. Later, the Dirichlet or the Neumann 
boundary conditions will be prescribed on T D and T N , respectively. Considering T > 0 a fixed 
time horizon, we set 


/ := [0, T], Q := (0,T)xl2, S := JxT, S D :=/xr D , S N :=/xr N . 


Further, and will denote the set of symmetric or symmetric trace-free (= deviatoric) 
(<ix<i)-matrices, respectively. For readers’ convenience, let us summarize the basic notation 
used in what follows: 


d = 2,3 dimension of the problem, 

Rf x ^ := {A <E R dxd ; A = A T }, 

:= {A e tr A = 0}, 

u : Q —> displacement, 

7T : Q —» R)jg V d plastic strain, 

C : Q — > [0,1] damage variable, 
a > 0 activation energy for damage, 

/:£ n —> R rf applied traction force, 
g:Q —> R d applied bulk force (as gravity), 
<t 0 i : Q —>■ Rfym elastic stress, 


e e i : Q —> Rf-jfm elastic strain, 
e = e{u) = e e i+vr = iVu T + total 

small-strain tensor, 
C : [0, lj -» R 34 elasticity tensor (dependent on £), 
H € R 3 hardening tensor (independent of C) 

S C Rjg V d the elastic domain (convex, hits' 3 0), 
w D : E d — > R^ prescribed boundary displacement, 
K\ > 0 scale coefficient of the gradient of plasticity, 
> 0 scale coefficient of the gradient of damage, 
b > 0 activation energy for possible healing. 


Table 1. Summary of the basic notation used through the paper. 


The state is formed by the triple q := (u, i r, £). Considering still a (small but fixed) regular¬ 
izing parameter £ > 0, the governing equation/inclusions read as: 

divcr e i + g = 0 with cr el = C(£)e e i and e e \ = e(u)— n, (momentum equilibrium) 
dd* s (Tr) 3 dev <T e i — Hbr + kiAtt, (plastic flow rule) 

^-o,6](C) ^ -^ c '(C)e e i : e el + k 2 div(| VCr _2 VC) - 2V [0il] (C), (damage flow rule) 
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(2.1a) 

(2.1b) 

(2.1c) 



with Ss the indicator function to S and its convex conjugate and with “dev” denoting the 
deviatoric part of a tensor, i.e. dev A := A — tvA/d. Here, [C(C)e]jj means Ylt i =1 Qjw(C) e w- 
We employed the regularizing term with a regularizing parameter £ > 0 with an exponent 
to be assumed suitably big, namely r > d. This regularization will facilitate analytical well- 
posedness of the problem and, because the gradient-damage term degenerates at VC = 0, 
its influence is presumably small if £ is small and VC not too large. 

Of course, (12. ip is to be completed by appropriate boundary conditions, e.g. 


u = w D on T D , (2.2a) 

Oei -n = f on T N , (2.2b) 

Vnn = 0 and VC -ft = 0 on T (2.2c) 

with n denoting the unit outward normal to O. We will consider an initial-value problem 
for (12.11) - fj2.2[) by asking for 

w(0) = m 0 , vr(0) = vr 0 , and C(0) = Co- (2.3) 


In fact, as u does not occur in (12. 1 p . m 0 is rather formal and its only qualification is to make 
S (0, Mo, tto, Co) finite not to degrade the energy balance (13.ldf) on [0, t]. 

After considering an extension u D = m d (t) of w D (t) from (I2.2al) on the whole domain H, 
it is convenient to make a substitution of u + m d instead of u into (12.ip — (I2.2p . we arrive to 
the problem with time-constant (even homogeneous) Dirichlet boundary conditions. More 
specifically, 


e e i in (12.lbh replaces by e e i = c(mTm d ) —7r, and (2.4a) 

w-o in (I2.2ap replaces by 0. (2.4b) 

Assuming (C(C)e(u D )-n = 0 on T N for any admissible (, this transformation will keep / in 
(I2.2bl) unchanged. 

Actually, (12. lbl) represents rather the thermodynamical-force balance governing damage 
evolution while the corresponding flow rule is written rather in the (equivalent) form 

7r G Ns(dev <Tei — Hbr — KiAn) with a e \ = C(C)e e i (2.5) 

and with Ns denoting the set-valued normal-cone mapping to the convex set S. An analogous 
remark applies to (I2.1cp . The system (12.ip with the boundary conditions (12.2p has, in its 
weak formulation, the structure of an abstract Biot-type equation (or here rather inclusion, 
cf. also e.g. HEBEI]): 

d^{q) + d q S(t, q) 9 0 (2.6) 

with suitable time-dependent stored-energy functional S and the state-dependent (pseudo)- 
potential of dissipative forces AS. Equally, as already used in (12.5p . one can write (12.6p as a 
generalized gradient flow 


q e d^AS* ( — dS{t, q)) 


(2.7) 
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where £%* denotes the conjugate functional. The governing functionals corresponding to 
m-m after the transformation (\2A\ are: 

7 T,C) := f ^C(C){e(u+u 0 (t))-Tr) : (e(u+it D (t))- n) + \m.tt : tt 

+ ^"|Vvr| 2 + —|VC| r + <5 [0 ,i](C) — g(t)-udx — / /(t)-udS, (2.8a) 

2 r ^r N 

^(7r, C) = &i{n) + ^2(C) := [ + a(~ + b( + dx (2.8b) 

Jo. 

where z + := max(z, 0) and := max(— z, 0) > 0. Note that the damage does not affect the 
hardening, which reflects the idea that, on the microscopical level, damage in the material 
that underwent hardening develops by evolving microcracks and even a completely damaged 
material consists of micro-pieces that bear the hardening energy ^Ebr : tt stored before. This 
model preserves coercivity of hardening even under a complete damage but the analysis 
below admits only incomplete damage. If ( 1 —y C(£)e:e is strictly convex for any e 7 ^ 0, we 
speak about a cohesive damage which exhibits a certain hardening effect so that the needed 
driving force increases when damage is to be accomplished. We can thus model quite a 
realistic response to various loading experiments, as schematically shown on Figured] for the 
case of a possible complete damage (whose analysis remains open, however). Note that, due 
to the “incompressibility” constraint tr tt = 0 , no plastihcation is triggered under a pure 
tension or compression loading. 



Figure 1: Schematic response of the stress cr el to the total strain e during a “one-dimensional” 
tension (left) or shear (right) loading experiment under a stress-driven scenario. The latter 
option combines plasticity with eventual complete damage. Dashed lines outline a response 
to unloading, C = C(() refers to Young’s modulus (left) or the shear modulus (right). 

Let us further note that [u, tt) H > <f>(t, u, n, () is smooth so that d q <f> = {^} x {(T)} x 
with S' u and S’) denoting the respective partial Gateaux derivatives and (12.6[) can thus be 
written more specifically as the system: 


<(t,M,7r,C) = 0, (2.9a) 

d^i(fr) + <^(f, u,n,C) 3 0, (2.9b) 

d& 2 (() + d^(t,u,TT,C) 9 0. (2.9c) 

Remark 2.1 (Irreversible damage in engineering models). Usual engineering models con¬ 
sider b = 00 , i.e. no healing is allowed. In fact, due to an essentially missing driving force 
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for healing, our modification b < oo would not have any influence on the evolution if it 
were not any VC-term in the stored energy. Thus, if the healing threshold b is big and the 
gradient-term coefficient k 2 > 0 is small, we expect to have essentially the (usually desired) 
unidirectional evolution as far as the damage concerns. 

Remark 2.2 ( Surface variant of the damage/plasticity). A similar scenario distinguishing 
tension (which leads to damage without plastification) and shear (with plastifying the mate¬ 
rial before damage) as in Figure |T] was used in a surface variant to model an adhesive contact 
distinguishing delamination in the opening and in the shearing modes, devised in [32[ [33] 
and later implemented by the fractional-step discretisation with checking the approximate 
maximum-dissipation principle in [25], 35,32]. An additional analytical difference is that, in 
contrast to our bulk model here, the surface variant allows for irreversible damage that does 
not need any gradient. 

Remark 2.3 (Other material models). A separately convex stored energy S(t, •) occurs also 
in other models. E.g., some phenomenological models for phase transformations in (poly¬ 
crystalline) shape-memory materials [3?j gives f the meaning of a volume fraction (instead 
of damage) and n a transformation strain (or a combination of the plastic and the trans¬ 
formation strains), and the total strain decomposes as e(u) = e e \ + (tt rather than (I2.1al) 
or makes n dependent on ( (which is then vector-valued). Considering the degree-1 homo¬ 
geneous dissipation potential, most of the considerations in this paper can be applied to 
such a model, too; in fact, the only difference would be the nonsmoothness of S also with 
respect to n variable. A similar (in general non-convex) model have been also considered 
in |3 im UM [22] EH] although sometimes special choices of elastic moduli leading to convex 
were particularly under focus while the dissipation is made state-dependent. 

3 LOCAL SOLUTIONS 

We will use the standard notation kF 1,p (r2) for the Sobolev space of functions having the gra¬ 
dient in the Lebesgue space L p (h2; M d ). If valued in M n with n > 2, we will write W l,p {fl\ M n ), 
and furthermore we use the shorthand notation R 1 (h2;M n ) = M n ). We also use the 

notation of “ • ” and “ : ” for a scalar product of vectors and 2nd-order tensors, respectively, 
and later also “ for 3rd-order tensors. For a Banach space A", L P (J; A) will denote the 
Bochner space of A-valued Bochner measurable functions «:/—>■ A with its norm ||u(-)|| 
in L P (J), here || • || stands for the norm in A. Further, W 1,P (I; X) denotes the Banach 
space of mappings u : I X whose distributional time derivative is in L P (J; A), while 
BV(/;A) will denote the space of mappings u : I X with a bounded variations, i.e. 
sup 0<to<tl< <tn _ 1<tn<T Y^i=i \\ u {li)~ u ifi- i)l| < °° where the supremum is taken over all fi¬ 
nite partitions of the interval / = [0, T\. By B(/;A) we denote the space of bounded 
measurable (everywhere defined) mappings / —y A. 

The concept of local solutions has been introduced for a special crack problem in [40] and 
independently also in [39], and further generally investigated in [2D]. Here, we additionally 
combine it with the concept of semi-stability as invented in [29]. We adapt the general 
definition directly to our specific problem, which will lead to two semi-stability conditions 
for ( and 7r, respectively: 

Definition 3.1 (Local solutions). We call a measurable mapping ( u , 7r, () : I —> R 1 (h2; R d ) x 
M^ v d ) x kF 1,r (r2) a local solution to the elasto-plasto-damage problem (12. Ill - (12.3|) if the 
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initial conditions (12.3j) are satisfied and, for some J C I at most countable (containing time 
instances where the solution may possibly jump), it holds that: 


Vte/\J: = 0, (3.1a) 

Vte/\J Wetf^njM^ev) : <?(*,«(*), 7r(t),C(*)) < <^(t, u(t), tt, C(*)) + ^i(7r-7r(£)), (3.1b) 

Vte/\{0} VCeW /1 ’ r (fi), 0<C<1: 

^(f,u(f),7r(f),C(t)) < £{t,u(t), 7r(£),C) +^’ 2 (C-CW), (3-lc) 

V0<ti<f 2 <T : <f’(t 2 ,M(t 2 ),7r(t 2 ),C(t2)) +Diss^ 1 (7r; [ii,i 2 ]) + Diss^ 2 (C; [ti,£ 2 ]) 

< ^(ti,it(ti),7r(ti),C(ti)) + f (3.Id) 


where 

Diss^(7r; [ r,s ]) := 


sup 

TVeN 


r<t Q < ''^T 1 ® 


* /■ 
j=i 


- 7 r(i,)) dx, and similarly (3.2a) 


iV 


Diss .« 2 (C; M) : = 


sup 

jven 


^ / a(7r(ij_i)-7r(ij)) + 6(vr(t i _ 1 )-vr(t i )) + dx. 




3 = 1 


(3.2b) 


Let us comment the above definition briefly. Obviously, (12.1 al) after transforming the 
boundary condition (\2A\i means precisely (13.lap , which more in detail here means that 
J n C(((t))(e(u(t)—w 0 (t)) — 7r):e(u)dx = J Q g-vdx + J^f-vdS for all v G H l { fl; M d ) with 
u|r D = 0, i.e. the weakly formulated Euler-Lagrange equation for displacement. Note that 
(13. lap specifies also the boundary conditions for u, namely u = 0 on r D because otherwise 
<f(t, u, 7r, () = 00 would violate (13. lap for v which satisfies v — 0 on r D , and also a ei ■ n = f 
on r N can be proved by standard arguments based on Green’s theorem. Equivalently, one 
can merge (13.lap with (13.lb[) to a single condition 

VteI\J \/(u,n)eH\n;R d ) x H\fl-,R d d ^), u\ Td = 0 : 

<f(t,ii(t),7r(t),C(t)) < &(t,u,n,C(t)) +&i(n-n(t))-, (3.3) 


which reveals that Definition 13. II just copies the concept of local solutions from J39| [JO] here 
generalized for the case of non-vanishing dissipation 7 ^ 0 . As is homogeneous degree- 
1, always C d&i(0) and thus (12.9bp implies d&i(0) + d^(u, 7r, £) 3 0. From the 

convexity of when taking into account that M\{0) = 0 , the latter inclusion is equivalent 

to ® x (v) + (d 7T S’(u(t),7r(t),((t)),v) > 0 for any v e H 1 ( 0;M^ v d ). Substituting v = z — z(t ) 
and using the convexity of S'ft, u, (, •), we obtain the semi-stability (13.lbf) of tt at time t. 
Analogously, we obtain also (13.1 cl) from (12.9cl) ; note that we do not require its validity at 
t = 0 so that we do not need to qualify the initial conditions as far as any (semi)stability 
concerns. Eventually, (13.1 dp is the (im)balance of the mechanical energy. Note that, in view 
of (I2.8a,p . the last term in (14.6dp involves 


7T,C)= / C(C)(e(u D )) : (e(u+u D )- 7 r) - g(t)-udx - j f(t)-udS. 

Jn J r N 

This is equivalent (or, if ${t, ■, ■, •) is not smooth, slightly generalizes) the standard definition 
of the weak solution to the initial-boundary-value problem (j2.ip ~ fj2.3p . cf. [31], Prop. 2.3] for 
details. 
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To be more precise, the concept of local solutions as used in [20, 30] requires J only 
to have a zero Lebesgue measure and also (13.lc[) is valid only for a.a. t. On the other 
hand, conventional weak solutions allow even (13.ldj) holding only for a.a. t\ and O- Later, 
our approximation method will provide convergence to this slightly stronger local solutions, 
which motivates us to have tailored Definition 13.11 straight to our results. 

Actually, local solutions form essentially the largest reasonable class of solutions for rate- 
independent systems as (I2.ip ~ fl2.3p considered here. It includes the mentioned energetic 
solutions m [22], the vanishing-viscosity solutions, the balanced-viscosity (so-called BV) 
solutions, parametrized solutions, etc.; cf. [201 2Tj for a survey, and also stress-driven-like so¬ 
lutions obeying maximum-dissipation principle in some sense, cf. Remark 13.21 The energetic 
solutions have often tendency to undergo damage unphysically early; cf. na for a compari¬ 
son on several computational experiments on a similar type of problem. The approximation 
method we will use in this article leads rather to the stress-driven option, cf. Remarks 13.21 
and 14.31 below. 

Remark 3.2 ( Maximum-dissipation principle). The degree-1 homogeneity of and <^2 
defined in (12.8b|) allows for further interpretation of the flow rules (12.9bp and (12.9cp . Using 
maximal-monotonicity of the subdifferential, (I2.9bj) means just that (£ — £ p i ast ,u — 7r) > 0 
for any v and any £ G (v) with the driving force £ p i ast = —<£%(t,u,n,Q■ Li particular, 
for v — 0, defining the convex “elastic domain” K\ = d&i(0), one obtains 

(Cpiast(t), n{t)) = max (£, n (£)) with some £piast (t) = u(t), vr(£), £(£)). (3.4a) 

CeiCi 

To derive it, we have used that £ p i ast € C d&i(0) = K\ thanks to the degree-0 

homogeneity of d&i, so that always (£ p i ast , ^) < max^ Ki (£ 7 f). The identity (I3.4al) says that 
the dissipation due to the driving force £ p i ast is maximal provided that the order-parameter 
rate 7r is kept fixed, while the vector of possible driving forces £ varies freely over all admissible 
driving force from K\. This just resembles the so-called Hill’s maximum-dissipation principle 
articulated just for plasticity in [16]. Also it says that the rates are orthogonal to the elastic 
domain Ad, known as an orthogonality principle |43j . Actually, R. Hill [TB] used it for 
a situation where <o{t, ■) is convex while, in a general nonconvex case as also here when 
damage is considered, it holds only along absolutely continuous paths (i.e. in stick or slip 
regimes) which are sufficiently regular in the sense 7r is valued not only in A 1 (f2;Mjg V d ) but 
also in H l (Q:R^ v d )* while it does certainly not need to hold during jumps. Analogously it 
holds also for £, defining K 2 := <9^(0), that 

(fdam(*),C(*)) = max (£,£(£)) with some £ dam (L) e -d^(t, u(t), n(t), ((t)). (3.4b) 

?6A' 2 

Here, d^(t,u,n,C) is set-valued and its elements should be understood as “available” 
driving forces not necessarily falling into K 2l while £d am £ K 2 is in a position of an “ac¬ 
tual” driving force realized during the actual evolution. As is smooth, the 

maximum-dissipation relation f!3.4al) written in the form (—£ > l(t,u(t),Tr(t),C(t)), / fr(t)) = 
max(A'i, 7r(t)) = ^i(7r(t)) summed with the semistability (13.lbj) which can be written in 
the form &i(ri) + (<sy(£, u(t), 7 r(f), £(£)), n) > 0 thanks to the convexity of $(t, u, •, £) yields 


^iW + {<?(*>CM ). 1 - i(t)) > 


( 3 . 5 ) 




























for any tt, which just means that £ p i ast (t) = —£ > f(t,u(t),n(t),C(t)) G d&i(n (f)), cf. f]2.9b[) . 
This exactly means that the evolution of n is governed by a thermodynamical driving force 
£piast (we say that it is “stress-driven”) and it reveals the role of the maximum-dissipation 
principle in combination with semistability. Using the convexity of <?(t, u, 7r, •), a similar 
argument can be applied for (13.4bl) in combination with semistability (13. left even if £{t, u, 7r, •) 
is not smooth. 


Remark 3.3 ( Integrated maximum-dissipation principle). Let us emphasize that, in general, 
7r and £ are measures possibly having singular parts concentrated at times when rupture 
occurs and the solution and also the driving forces need not be continuous. Even if d and £ 
are absolutely continuous, in our infinite-dimensional case the driving forces need not be in 
duality with them, as already mentioned in Remark 13.21 So (13.41) is analytically not justified 
in any sense. For this reason, an Integrated version of the Maximum-Dissipation Principle 
(IMDP) was devised in [31] for a bit simpler case involving only one maximum-dissipation 

relation. Realizing that rriax^^ (£, 7r) = and similarly max~ eff2 (^, £) = ^ 2 (C)> the 

integrated version of (13.4|) reads here as: 


r^2 


rt2 


/ fpiast(i) dvr(t) = / ^i(7i)d t 
Jti Jt\ 

AdamWdC(t)= [%(()dt 

hi Jt x 


with £piast(£) = -<^(t,u(t),7r(t),C(t)), (3.6a) 

with some £ da m(£) G-0 c <?(t, u(t), 7r(t), £(t)) (3.6b) 


to be valid for any 0 < t\ < t 2 < T. This definition is inevitably a bit technical and, without 
sliding too much into details, let us only mention that the left-hand-side integrals in (13.61) 
are the lower Riemann-Stieltjes integrals suitably generalized, defined by limit superior of 
lower Darboux sums, i.e. 


[ f(f) Mt) 

J r 


N 

limsup inf (£(£), 1))> 

n gn te[tj-i,tj] 

r=t 0 <t 1 <...<t N _ 1 <t N =s J 


(3.7) 


relying on that the values of £ are in duality with values of z (but not necessarily of z) and on 
that the collection of finite partitions of the interval [r, s] forms a directed set when ordered 
by inclusion so that “limsup” in (13.71) is well defined. Let us mention that the conventional 
definition uses “sup” instead of “limsup” but restricts only to scalar-valued ( and z with z 
non-decreasing. The limit-construction (13.7[) is called a (here lower ) Moore-Pollard-Stieltjes 
integral [231 [27 ] used here for vector-valued functions in duality, which is a very special case 
of a so-called multilinear Stieltjes integral. Like in the mentioned classical scalar situation of 
lower Riemann-Stieltjes integral using “sup” instead of “limsup”, the sub-additivity of the 
integral with respect to u and to v holds, as well as additivity with respect to the domain 
holds. 

The right-hand-side integrals in (13.61) are just the integrals of measures and equal to 
Diss^ (tt; [£i,£ 2 ]) and Diss^ 2 (£; [ti, t 2 ]), respectively. Equivalently, in view of the definition 
(13.2|) . they can be also written as d7r(f) and 2 d£(f), where the integrals can again 

be understood as the lower Moore-Pollard-Stieltjes integrals (or here even simpler as the 
mentioned lower Riemann-Stieltjes integrals) modified for the case that the time-dependent 
linear functionals £ are replaced by nonlinear but time-constant and 1-homogeneous convex 
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functionals M’s. Alternatively, though not equivalently, denoting the internal variables z = 
(•7T, £), the IMDP (13.6p can be written “more compactly” as 

H 2 rt2 

/ £(t) dz(t) = / & dz(t) with some £(t)E~d z £’(t,u(t),z(t)). (3.8) 

Jti Jt\ 

Both IMDP (13.61) or (13.81) are satished on any interval [ti,t 2 ] where the solution to 
(12.9|) is absolutely continuous with sufficiently regular time derivatives; then the integrals 
in (13.6p are the conventional Lebesgue integrals, in particular the left-hand sides in (13.61) 
are (£ p i ast (t ), 7r(t)) dt and (£dam(£), C(t)) di, respectively. The particular importance of 
IMDP is especially at jumps, i.e. at times when abrupt damage possibly happens. It is 
shown in I2UEH on various finite-dimensional examples of “damageable springs” that this 
IMDP can identify too early rupturing local solutions when the driving force is obviously 
unphysically low (which occurs quite typically in particular within the energetic solutions of 
systems governed by nonconvex potentials like here) and its satisfaction for left-continuous 
local solutions indicates that the evolution is stress driven, as explained in Remark 13.21 On 
the other hand, it does not need to be satisfied even in physically well justified stress-driven 
local solutions. For example, it happens if two springs with different fracture toughness 
organized in parallel rupture at the same time, cf. [21]) Example 4.3.40], although even in 
this situation our algorithm (14.21) below will give a correct approximate solution, cf. Figure [6] 
below. Therefore, even the IMDP (|3.6[) may serve only as a sufficient aposteriori condition 
whose satisfaction verifies the obtained local solution as a physically relevant in the sense 
that it is stress driven but its dissatisfaction does not mean anything. Eventually, let us 
realize that, as a consequence of the mentioned definitions, we have 


-*2 


rt2 


rt2 


fpiastfa) dvr(t) + / fdam(*)dC(t) < / f(t) d(vr, £) (t) and 


'tl 




rt2 


rt2 


rt 2 


'ti 


?idn (t) + / & 2 d( (t) = / (&i+& 2 ) d(7T, C)(t)- 


Iti 




(3.9a) 

(3.9b) 


As there is only inequality in (13.9a[) . the IMDP (13.81) is less selective than (13.6|) in general. 
Moreover, we will rely rather on some approximation of IMDP, cf. Remarks 14.31 and 16.21 
below. 


4 SEMI-IMPLICIT TIME DISCRETISATION AND ITS CON¬ 
VERGENCE 

To prove existence of local solutions, we use a constructive method relying on a suitable 
time discretisation and the weak compactness of level sets of the minimization problems 
arising at each time level. When further discretised in space, it will later in Sect. 0 yield a 
computer implementable efficient algorithm. Let us summarize the assumption on the data 
of the original continuous problem: 

C(-),He ^dxdxdxd positive definite, symmetric, C : [0,1] ->• R dxdxdxd continuous, (4.1a) 
a, b, >0, Sc convex, bounded, closed, lilts' 9 0, (4-lb) 

w D e VF 1,1 (/; W 1 / 2 ’ 2 (r D ; M d )), (4.1c) 

g ewCi;L^)) with p{l\ dl(d+2) £||>s (4-id) 
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with p {= l- 2 /d fold is! < 41e > 

K? To,Co) e (4.1f) 

The qualification (14.lc[) allows for an extension w D of w D which belongs to W 1}1 (I ; M d )); 

in what follows, we will consider some extension with this property. 

For the mentioned time discretisation, we use an equidistant partition of the time interval 
/ = [0, T ] with a time step r > 0, assuming T/t E N, and denote {u k }^ 0 an approximation 
of the desired values u(kr), and similarly 0 is to approximate ((kr), etc. 

We use a decoupled semi-implicit time discretisation with the fractional steps based on 
the splitting of the state variables governed by the separately-convex character of S’(t,- 
This will make the numerics considerably easier than any other splitting and simultaneously 
may lead to a physically relevant solutions governed rather by stresses (if the maximum- 
dissipation principle holds at least approximately in the sense of Remark 14.31 below) than 
by energies and will prevent too-early debonding, as already announced in Section [3j More 
specifically, exploiting the convexity of both <£(t, and S’ft, u, tt, •) and the additivity 

& = &i{k) T (C) > this splitting will be considered as (u , 7r) and (. This yields alternating 
convex minimization. Thus, for (tt* - 1 ,^ -1 ) give 11 ? we obtain two minimization problems 

minimize (u, 7T, Cr” 1 ) + ^\{tt— ir k ~ l ) 
subject to (m,7t) G H x M^* v d ), u\ Fo = 0, 

with S J f := ^{kr, •, •, •) and, denoting the unique solution as ( u k , tt*), 

minimize (ofivJt , ir k Q + ^(C - Cr _1 ) 

subject to (6 W 1>r (Q), 0 < C < 1, 


(4.2b) 


(4.2a) 


and denote its (possibly not unique) solution by (f. Existence of the discrete solutions 
{u k , 7 (f) is straightforward by the mentioned compactness arguments. 

We define the piecewise-constant. interpolants 

u T (t) = u k & u T (t) = it* -1 , 

TT T (t) = 7T k & 7 r_ T (t) = 7T* 1 , 

Ut) = C & C T (t) = Cf~\ 

&’r(t,U,n,C) = <%T 0?7bC) 

Later in Remark 14.31 we will also use the piecewise affine interpolants 


> for (k—l)r <t < kr. (4.3) 


TT T (t) = + ^7T k r -\ 

Cr{t) = 


for (k—l)r < t < kr. 


(4.4) 


The important attribute of the discretisation (j4.2[) is also its numerical stability and 
satisfaction of a suitable discrete analog of (13. ip . namely: 

Proposition 4.1 (Stability of the time discretisation). Let (14. ip hold and, in terms of 
the interpolants (14. 3p . (u T ,7t T ,( T ) be an approximate solution obtained by (14. 2p . Then, the 
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following a-priori estimates hold 

II u " 


< (1 

II7F II < a 

lr lr llL°°(7;_ff 1 (Q;R^ v li ))nBV(7;L 1 (r2;R^ v d )) — ’ 

ll^ T IL°°(n)nBV(/;L 1 (n)) — C- 


(4.5a) 

(4.5b) 

(4.5c) 


Moreover, the obtained approximate solution satisfies for any t G 7\{0} the (weakly formu¬ 
lated) Euler-Lagrange equation for the displacement: 

£((t T ,u T (t),T? r (t),( T (t)) =0, (4.6a) 

with t T \— min{/cr>f; k&N}, two separate semi-stability conditions for ( T and tt t : 

WeH 1 (£l)W$£) : £(t T ,u T (t),n T (t),( T (t)) < #(t T ,u T (t),7f,( T (t)) +^i(5f-7r T (t)), (4.6b) 

vcw 1 >r (fi), o<C<i: 

£{t r , U r (t),W T (t),( T (t)) < S(t T ,U T {t), 7f r (t),C) + & 2 {(~Ct( t)), (4.6c) 

and, for all 0 < ti < t 2 < T of the form ti = kiT for some ki gN, the energy (im)balance: 

^(t2,u T (t2),7r T (t2)Xr(h)) + Diss^ (tt t ; [ti,t 2 ]) + Diss^ 2 (C T ; [* 1 ,^ 2 ]) 

< + f ^(t,U T (t),7f(t)X(t))dt. (4.6d) 

Jtx 

Sketch of the proof. Writing optimality condition for (I4.2a| in terms of u, one arrives at 
(I4.6al) . and comparing the value of (j4.2aj) at (u k r ,Tr , f) with its value at (u 1 ), tt) and using the 
degree-1 homogeneity of one arrives at (I4.6bl) . 

Comparing the value of (14.2b[) at Cf) with its value at ( and using the degree-1 homogeneity 
of & 2 , one arrives at (I4.6cj) . 

In obtaining ( I4.6djl . we compare the value of (14. 2 all at the minimizer (u 1 ), tt!)) with the 
value at (u^T 1 , t^ 1 ) 1 ) an d R ie va l ue °f H4.2b[) at the minimizer (f with the value at _1 and 
we benefit from the cancellation of the terms ±S , (kT,u l (,TT l ),(f~ 1 ). We also use the discrete 
by-part integration (= summation) for the <f/-term. 

Then, using (14.6dj) for ti = 0 and the coercivity of S’it, ■, •) due to the assumptions 

m, we obtain also the a-priori estimates (14.51) . □ 

The cancellation effect mentioned in the above proof is typical in fractional-step methods, 
cf. e.g. [201 Remark 8.25]. Further, note that (I4.6j) is of a similar form as (13. 1 p and is thus 
prepared to make a limit passage for r —* 0: 

Proposition 4.2 (Convergence towards local solutions). Let (14.ip hold and let (u T , 7 f T , £ r ) be 
an approximate solution obtained by the semi-implicit formula (14.211 . Then there exists a sub¬ 
sequence (indexed again by r for notational simplicity) and u G B([0, T\; i7 1 (0; R d )) and tt G 
B([0,T]; Wffl; Rto))rtBV([0,71; R^)) and( e B([0,T];H' 1 -’'(SJ))nBV([0,T]; L'(Q)) 

such that 


u T (t) —> u(t) 

TT r {t) — > TT(t) 

C r(t) ~> C (t) 


mH 1 (n-,R d ) for all t G [0, T\, 

in H 1 (LI; M^ v d ) for all t G [0, T\, 

in W l,r (Lf) for all t G [0, T}. 


(4.7a) 

(4.7b) 

(4.7c) 
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Moreover, any (u, 7r, () obtained by this way is a local solution to the damage/plasticity 
problem in that sense of Definition ^. 1\ 


Proof. By a (generalized) Helly’s selection principle, cf. also e.g. [201 2T], we choose a subse¬ 
quence and 7T G B([0,T];i/ 1 (^^ e x v d ))nBV([0,r];L 1 (fi;M^))andC, C* e B([0, T]; VF 1,r (Q)) 
n BV([0,T];L 1 (fi)) so that 

7f r (f) —^ 7r(f) in Mj X v d ) for all fe[0,T], (4.8a) 

Cr(t)^({t) & C T (t)^(*{t) in W 1,r (Pt) for all te[0,T\. (4.8b) 

Now, for a fixed t E [0, T], by Banach’s selection principle, we select (for a moment) further 
subsequence so that 


u T (t) —*■ u(t) in H 1 (Q;R d ). (4.9) 

We further use that u T {t) minimizes S’(t T , •, 7f T (t), £ (£)) with t T := niin{/cr > f; k E N}. 
Obviously, t T —> t for r —> 0 and, by the weak-lower-semicontinuity argument, we can easily 
see that u(t) minimizes the strictly convex functional ${t, •, £*(£), 7t(£)); this is indeed simple 
to prove due to the compactness in both n and ( due to the gradient theories involved. Thus 
u(t) is determined uniquely so that, in fact, we did not need to make further selection of a 
subsequence, and this procedure can be performed for any t by using the same subsequence 
already selected for (14. 8p . Also, u : [0, T] —> H l [ Q;R d ) is measurable because 7r and £* are 
measurable, and S l f l (t,u(t),'ir(t),C*(t)) — 0 for all t. 

The key ingredient is improvement of the weak convergence (14.81) and (I4.9[) for the strong 
convergence. For the strong convergence in u and 7r, we use the uniform convexity of the 
quadratic form induced by C(£), H, and K\ with the information we have at disposal from 
(14.6bj) leading, when using the abbreviation e el = e{u—u D ) — 7r and e elr = e(u T —u DtT ) — n T , 
to the estimate: 

[ C(C T (t))(e ehT (t)-e el (t)) : (e AiT (t)-e A (t)) 

Jn 

+ El(7f r (t)—7r(t)) : (jf T (t) 7r(t)) + ^ | VT r (t)-V7r(t) | 2 dx 

< /-C(C (t))e e i (t) : (e eljT (t)—e el (t)) - (Htt( t)-| T (t)) : (7r T (t)-7r(t)) 

Jn 

+ y V7r(f) : v ( 7 f T (t)-n(t)) - J T (t)-(u T (t)-u(t)) d x - f g T (t)-(u T (t)-u(t)) dS -> 0 
z 4r N 

where we use some f T (t) E dSg(W T (t)) which solves at time t in the weak sense the discrete 
plastic flow-rule + Hl7f T — devoy = k,iAtc t with a T = C(£ )e eliT . Thus we proved 

ejt) strongly in L 2 (M,R d //) 

together with (]4.7b|l . Realizing that e(u T (t )) = e(tZ DiT (f)) + 7f T (t) + e elr (t), we obtain also 
e(u T (t)) —> e(u(t )) strongly in L 2 (M,R d //), and thus also (I4.7al) . Note that we exploited 
the gradient theory for plasticity which ensures that the sequence (^ T ) r >o; which is bounded 
in L°°(Q; R d d f v d ) because the plastic domain S C R){ x v d is bounded, is relatively compact in 
so that the term ■ (7r r (t)—7r(t)) dx indeed converges to zero because 

7T T (0^7T(t) in H\n- R^). 
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The convergence (j4.7cj) can be proved by using the uniform-like monotonicity of the set¬ 
valued mapping ( i-^ cM[ 01 ](C) — k 2 div(| VC| r_2 VC) : W 1,r (Q) =1 IU 1,r (r2)*. Analogously to 
(12. left , we can write the discrete damage flow rule after the shift (12.4ft as 

£dam,r + : e d>T = K 2 div(I'VCr| r_2 'V(r) “ Vr (4.10a) 

with some £ dam , r e <9<5*_ q6] (C r ) and rj T edS [01] (( T ) (4.10b) 

with the boundary condition V£ r ■ n — 0 on E; in (14.101) . Cdam,r and rj T are considered 
piece-wise constant in time, consistently with our bar-notation. An important fact is that 
Cdam , T (t) is valued in [—6, a] and hence a-priori bounded in L°°(O); here we vitally exploited 
the concept of possible (small) healing allowed. We can rely on Cdam,r(^) —*■ Cdam(^) in L°°(fl) 
for some ^-dependent subsequence and some Cdam(^)- Using that C'(C (^))e elT (t) : e el (i) is 
bounded and, due to (|4.7h .b). even has been proved converging in L 1 (h2) which is a subspace 
of VF 1,r (Q)* because r > d is considered. By the standard theory for monotone variational 
inequalities, we can pass to the limit in (14. lOj) at time t to obtain, in the weak formulation, 

£dam(f) + C'((,(t))e, I (t):e eI (t) = n 2 div(|VC(i)| r ~ 2 VC(i)) - 1 (<) with 7j(t)e85 M (C(t)). 

(4.11) 


Then, at any t, we can estimate 

k 2 limsup (IIVCtWII^LkcI)- l|VC(t)||^ ( 1 n . Rd) )(||VCr(t)|| L r (n;R;i) - ||VC(t)|| L r (n . Rd) ) 

k—yoo 

< iimsup f ^(ivcddr-w^w-ivcwr-wcwj'V&w-cw) 

k —>oo Jil 

= lim f C'(C T (i))e 8l , T (t) : e 8liT (i)(<r«-C«) - «2 |VC(t)r 2 VC(t)-V(C r (t)-C(t)) 

k— yoo J q 

- (Cdam (t)+v(t))(CAt)-C(t)) dx = 0 (4.12) 


where the last equality has exploited (14.lip . The important fact used for (14.121) is that 

c, (C T ( t ))®d,r(*) : eei,r(^)(Cr(t)-C(^)) 0 weakly in T 1 (0); (4.13) 

in fact, this convergence is even strong when realizing that £ T (t) —>■ C{t) i n L°°(Q), for 
which again r > d is exploited. From this, (j4.7c[) follows. Thus, from (14.12[) we can see 
that || VCr(^)||Lw(Q. R d) —y ||VC(0IU r (n ; K d ) an d, from uniform convexity of the Lebesgue space 
L r (Q]R d ), we eventually obtain (I4.7cl) . Actually, the specific value £dam (t) of the limit of (a 
f-dependent subsequence of) {^dam,r(^)}r>o which is surely precompact in IU 1,r (fl)* is not 
important and thus (14.7c() holds for the originally selected subsequence, too. 

Having the strong convergences (14.7ft proved, the limit passage from (14.6ft towards (13.11) 
is simple. In particular, by continuity of both BV-functions ((■) and £*(') on [0, T]\J for 
some at most countable set J, we have also £*(£) = ((t) at any t except at most countable 
the set J. □ 


Remark 4.3 ( Approximate maximum-dissipation principle). One can devise the discrete 
analog of the integrated maximum-dissipation principle (13.61) straightforwardly for the left- 
continuous interpolants (14.31) . required however to hold only asymptotically. More specifi¬ 
cally, in analog to (13.6ft formulated equivalently for all [0, i] instead of [ti,t 2 ], one can expect 
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an Approximate Maximum-Dissipation Principle (AMDP) in the form 

_ ? _ _ , _ _ 

CplastjT d.7T r ~ DiSS^^^; [0,t]) with £plast,r = ~ |_^rj 7r (-,M r ,7T r ,C T ), (4.14a) 

l Cdam,r d(^ r ~ DlSS^ 2 ((^ r , [0,t]) for SOIIie ()dam,T £ S T ( ', U T , 7T t j Ct) ; (4.14b) 

or, analogously to (13. 8 ft , 

£ T dz T ~ Diss^(2 r ; [0,t]) with some £ r G { — [ < ^-] 7r (-, w T , 7f r , C r )} x -d ( ~S‘ T (-,u T , n T ,(r), 

(4.15) 

where the integrals are again the lower Moore-Pollard-Stieltjes integrals as in (13.61) and where 

S’ri', u, vr, () is the left-continuous piecewise-constant interpolant of the values S{kr, u, 7r, £), 

? 

k = 0,1,..., T/r. Moreover, in (14.141) means that the equality holds possibly only asymp¬ 
totically for r —y 0 but even this is rather only desirable and not always valid. Anyhow, 
loadings which, under given geometry of the specimen, lead to rate-independent slides where 
the solution is absolutely continuous will always comply with AMDP (I4.14p . Also, some 
finite-dimensional examples of “damageable springs” in [2Tj [3Tj show that this AMDP can 
detect too early rupturing local solutions (in particular the energetic ones) while it generi- 
cally holds for solutions obtained by the algorithm (14.21) . Generally speaking, (14.14ft should 
rather be a-posteriori checked to justify the (otherwise not physically based) simple and nu¬ 
merically efficient fractional-step-type semi-implicit algorithm (I4.2[) from the perspective of 
the stress-driven solutions in particular situations and possibly to provide a valuable infor¬ 
mation that can be exploited to adapt time or space discretisation towards better accuracy in 
(14.14ft and thus close towards the stress-driven scenario. Actually, for the piecewise-constant 
interpolants, we can simply evaluate the integrals explicitly, so that AMDP (14.15ft reads 




K 


k =1 


v + „( Cr ‘- C y i )-+ b( c‘-c T ‘-') + di ? 

<£t\0 (4.16) 

where = -[^(uJ^.Cf- 1 ) and ^ e - 9 ^( 4 ,^,^) 


for some e T \ 0, where K = maxjfceN; kr < t}. Notably, in contrast to (13. 6 p and (13.8p . 

the AMDP (14.14ft and (14.151) are equivalent to each other as the limsup’s (cf. the definition 

(13.7ft ) in all involved integrals is attained on the equidistant partitions with the time step 

r and the “inf” in the Darboux sums is redundant. Evaluating the dualities, (I4.16P can be 

? 

written more explicitly as J n dx < e T \ 0 with the residuum 


K 


R* := 


k =1 ^ 

cK^Ky-'-ety-^^- 1 )) + Ebi-p 1 ): (y-y- 1 ) 

lc'(cf- 1 )(e(«;-'+ u ‘y)- ) r;-‘) ; (e(^- i +<y)- J r;- i )+e^,d(c;-c‘- 1 ) 


- KrVwp 1 : vpj-y- 1 ) - K 2 |vcPT VcP'-v^-cP 1 )) ( 4 . 17 ) 
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with some multiplier ^ onst T e ^[ 01 ](Cr) and with 0 2 for k = 1 equal to £ 0 - Note that 
Rfi cannot be guaranteed non-negative pointwise on 0, only their integrals over hi are non¬ 
negative. One can a-posteriori check the residua depending on t or possibly also on space, 
cf. also [351 142] for a surface variant of such a model or Figures [4H7] below. 


5 IMPLEMENTATION OF THE DISCRETE MODEL 


To implement the model computationally, we need to make a spatial discretisation of the vari¬ 
ables from the semi-implicit time discretization of Section 4. Essentially, we apply conformal 
Galerkin (or also called Ritz) method to the minimization problems (j4.2a.|) and (I4.2bl) which 
are then restricted to the corresponding finite-dimensional subspaces. These subspaces are 
constructed by the finite-element method (FEM), and the solution thus obtained is denoted 
by 

_ ( n ,k _k /-k \ 

Qrh ■— \ U rhi n rhiS,Th)i 

with h > 0 denoting the mesh size of the triangulation, let us denote it by 5^, of the domain 
0 considered polyhedral here. By this way, we obtain also the piecewise constant and the 
piecewise affine interpolants in time, denoted respectively by u T h and u T h , tt T h and 7r Th , and 
eventually Q r h and Crh- The simplest option is to consider the lowest-order conformal FEM, 
i.e. Pl-clcments for u, (, and ^ r. In Sect. [HI only the case d — 2 will be treated, so the 
previous analytical part have required r > 2 and we make an (indeed small) shortcut by 
considering r = 2. Moreover, we will not consider the loading on T N so we put / = 0. 

The material is assumed isotropic with properties linearly dependent on damage. The 
isotropic elasticity tensor is assumed as 


C ijki(C ) [(Ai—Ao)C + Ao]M« T [(/^i —/ x o)C + ho] {fiik^ji J rdii8jk) (5-1) 

where Ai,/ii and A 0 ,ho are two sets of Lame parameters satisfying Ai > A 0 > 0 and ji\ > 
/io > 0. Here, 6 denotes the Kronecker symbol. This choice implies that the elastic-moduli 
tensor is positive-definite-valued (and therefore invertible). The elastic domain S is assumed 
to satisfy 

s = W<<r Y }. ( 6 . 2 ) 

where <r Y > 0 is a given plastic yield stress. More specifically, the minimization problems 
(14. 2 j 1 after spatial discretisation rewrite as 


(u 


rhi ' rh 


7OJ = argmin 

u£W 1,o ° (n ; R d ) Jn 

dxc 
dev 


0 C (Cfh 1 )(e(M+t4,rh)- 7r ) : (e( M + M D,rh)-7r) 




«1 


u,tt elementwise affine on "h 2^^T 2 d-rh ^ 1 "b ^yI^" ^- 


fc-1 
rh I 


C rh = argmin / 

Cew 1 - 00 ^) Ji 

0<C<1on O 

(' elementwise affine on ^ 


Q^(0 { e ( U Th+ U D,Th)—' K Th) : {e{Urh+Uu,rh)-^Th) 


+ f|VCI 2 + a(c-cL 1 )- + KC-CSrT 


da;, (5.3a) 


da;. (5.3b) 


The damage problem (15,3bh represents a minimization of a nonsmooth but strictly convex 
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functional. To facilitate its numerical solution, we still modify it a bit, namely 


argmm 

C,C A ,C v ew 1 ’ oo (0) 
0<C<1 on n 

C,C A ,C V elementwise affine on 


n 


;C(()(e(u* h +u^ Th )-n* h ) : 


+ f|VC | 2 + <+f>C']dl, 


where ( A = (C — C r h 1 ) + and C V = (C — C th *) at all nodal points. 


(5.4a) 

(5.4b) 


We used additional auxiliary ‘update’ variables ( A and £ v which are also considered as 
Pl-functions. This modification can also be understood as a certain specific numerical 
integration applied to the original minimization problem (15.3b|) . It should be noted that 
( and Cr/7 1 are Pl-functions and, if we would require (I5.4bh valid everywhere on O, ( A 
and <C V could not be Pl-functions in general on elements where nodal values of (—• 
alternate signs. The important advantage of f!5.4bl) required only at nodal points while at 
remaining points it is fulfilled only approximately (depending on h ) is that fl5.4a,p actually 
represents a conventional quadratic-programming problem (QP) involving the linear and the 
box constraints 


C = (rh 1 


+ ( A ~C 


o < c A < i 


/•k —1 
rh J 


o < C v < C 


rh 


(5.5) 


A convex quadratic cost functional of this QP problem has only a positive-semidehnite 
Jacobian, since there are no Dirichlet boundary conditions on the damage variable (. Note 
that the optimal pair (£ A , ( v ) must satisfy () A (( V =0 in all nodes, i.e. both variables cannot be 
positive. This can be easily seen by contradiction: If () A <C V >0 in some node, then a different 
pair (C A — min{<C A , £ v }> C v — min{C A , C v }) would again satisfy the constraints (15.5|) but would 
provide a smaller energy value in fl5.4ap . 

As we have a-priori bounds of £ in hP 1,r (f2) uniformly in t, r, and h also if the modified 
problem (I5.4bp is considered (disregarding that we used r = 2 above), we have estimates 
also in Holder spaces also for ( A and and can show that the constraints (I5.4bp are valid 
everywhere on in the limit for h —> 0. Thus, an analogy of Proposition 14.21 for a successive 
limit passage h —> 0 and then r —> 0 might be obtained, although it does not have much 
practical importance for situations when ( h , r) (0, 0) simultaneously. 

A similar modification can be used also for fj5.3ap . In addition, one can then exploit the 
structure of the cost functional being the sum of a quadratic functional and a nonsmooth 
convex functional with the epi-graph having a “ice-cream-cone” shape. After introduction 
of auxiliary variables at each element, it can be transformed to a so-called second-order cone 
programming problem (SOCP), cf. [2TL Sect. 3.6.3], for which efficient codes exist. 

Other way is to use simply the quasi-Newton iterative method. This option was used 
also here. 


6 ILLUSTRATIVE 2-DIMENSIONAL EXAMPLES 

Finally, we demonstrate both the relevance of the model together with the solution concept 
from Sect .[2}{3] and the efficiency and convergence of the discretisation scheme from Sect. [4] 
together with the implementation from Sect. 0 on a two-dimensional example. 

The material: We consider an isotropic homogeneous material with the elastic properties 
given by Young’s modulus E Youn ~ 27GPa and Poisson’s ratio v = 0.2 in the non-damaged 
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state, which means that the elastic-moduli tensor in the form (15. ip takes Ai = 7.5 GPa 
and /i i = 11.25 GPa, while the damaged material uses lO'-times smaller moduli, i.e. Ao = 
750Pa and go = 112.5Pa in (15.ip . The yield stress from (15.2p and the kinematic hardening 
parameter are chosen as a Y = 2 MPa and HI = (i7 Young /20)I. The activation energy for 
damage is a = 1.2 kPa and the damage length-scale coefficient is k 2 = 0.001 J/m; the healing 
(used before for analytical reasons) was effectively not considered, cf. Remark 12.11 




Figure 2: Geometry of a 2-dimensional square-shaped specimens subjected to two tension¬ 
loading experiments; the right-hand side of the rectangle 11 combines Dirichlet condition in 
the horizontal direction and homogeneous Neumann condition in the vertical direction. 

The specimen and its loadings: We consider a 2-dimensional square-shaped specimen sub¬ 
jected to two slightly different loading regimes. Both of them consist in a pure “hard-devise” 
horizontal load by Dirichlet boundary conditions with the left-hand side T D fully fixed while 
the right-hand side r D /r N combines time-varying Dirichlet condition in the horizontal direc¬ 
tion with the Neumann condition in the vertical direction. The only (intentionally small) 
difference is in keeping a small bottom part of this vertical side free (see Fig.[23-left) or not 
(see Fig. [fright). As our model is fully rate-independent, the time scale is irrelevant and 
we thus consider a dimensional-less process time t E [0, 80] controlling the linearly growing 
hard-devise (= Dirichlet) load until the maximal horizontal shift 80 mm of the right-hand 

side r D /r N . 

The discretisation : In comparison with Section [5l we dare make a shortcut by neglecting the 
gradient term Vtt in the stored energy (I2.8aj) by putting K\ — 0, which allows for using only 
PO-elements for n. It also allows for transformation of the cost functional of (I5.3ap to a func¬ 
tional of the variable u only by substituting the elementwise dependency of 7r on u, see mm 
for more details. Then, the quasi-Newton iterative method mentioned in Section [5] is applied 
to solve u* h while 7is reconstructed from it. More details on this specific elasto-plasticity 
solver can be found e.g. in mmm- Here, the spatial P1/P0 FEM discretisation of the 
rectangular domain D uses a uniform triangular mesh with 2304 elements and 1201 nodes. 
The code was implemented in Matlab, being available for download and testing at Matlab 
Central as a package Continuum undergoing combined elasto-plasto-damage transformation, 
cf. [Hj. It is based on an original elastoplasticity code related to multi-threshold models [6], 
here simplified for a single-threshold case. It partially utilizes vectorization techniques of [28] 
and works reasonably fast also for finer triangular meshes. In contrast to the fixed spatial dis¬ 
cretisation, we consider three time discretisation to document the convergence (theoretically 
stated only for unspecified subsequences in Proposition 14.2p on particular computational ex¬ 
periments. More specifically, we used three time steps r = 1, 0.1, or 0.01, i.e. the equidistant 
partition of the time interval [0, 80] to 80, 800, or 8000 time steps, respectively. 
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Figure 3: Evolution of averaged elastic von Mises stresses , j dev cx e i| da; over time for the two 
experiments from Fig.\2\for three different time discretisations. The left experiment ruptures 
earlier under less stretch and leaving less plastihcation and remaining stress comparing to 
the right experiment. In both experiments, the convergence theoretically supported by 
Proposition 14.21 is well documented. 


Simulation results: The averaged stress/strain (or rather force/stretch) response is depicted 
in Figure 0 Notably, after damage is completed, some stress still remains (as is nearly 
independent on further stretch because the elastic moduli Ao and p 0 are considered very 
small). These remaining stresses are caused by non-uniform plastihcation of the specimen 
during the previous phases of the loading. One can also note that Figure 0-right imitates 
quite well the scenario from Figure [f|-right while Figure 01eft is rather a mixture of both 
regimes from Figure |T] and, interestingly, the rupture proceeds in three stages. The respective 
spacial distribution of the evolving state variable is depicted at few selected instants on 
Figures 0 and 0 It is well seen how a relatively small variation of geometry in Figure 0 
dramatically changes the spatial scenario and triggers damage in very different spots of the 
specimen. This is an expected notch-effect causing stress concentration and relatively early 
initiation of cracks at such spots, i.e. here such a notch is the point of the transition T N to 
r D /r N in Figure 0 left. The AMDP suggested in Remark 14.31 is depicted in Figures [6] and [71 It 
should be emphasized that the maximum-dissipation principle (as devised originally by Hill 
[IS]) is reliably satisfied only for convex stored energies as occurs during mere plastihcation 
phase, as also seen in Figure [3 while in general it does not need to be satisfied even in 
obviously physically relevant stress-driven evolutions, as already mentioned in Remark 13.31 
and which can be expected even here during massive fast rupture of a wider regions (but in 
spite of it, Figure 0 shows a good satisfaction of AMDP even during such rupture phases 
and in some sense demonstrate a good applicability of the model and solution concept and 
its algorithmic realization). 

Remark 6.1 (Symmetry issue). Actually, one could understood the square 1 x 1 in Fig. 0 
as one half of a rectangle with sides 2x1 with the right-hand side of the lxl square 
being the symmetry axis of the 2x1 rectangle which is then loaded from the vertical sides 
fully symmetrically. Engineers actually most routinely assume that such symmetry of this 
geometry would be inheritted by all (or at least by one) solution(s) and use the reduced 
geometries on Fig. 0 for calculations of the full 2 x 1-rectangle. We intentionally did not use 
this interpretation because, in fact, one can only say that the set of all solutions inherits the 
(possible) symmetry of the specimen and its loading but not particular solutions, and even 
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Figure 4: Evolution of spatial distribution of the state (u, n, () with also the von Mises stress 
and the residuum R from (j4.17[) at (equidistantly) selected instants for the asymmetric 
geometry from Fig.\2\-left. The deformation is visualized by a displacement u magnified 
250 x, and r — 0.1 was used. Damage occurs relatively early on the right-hand side due to 
the stress concentration and propagates in several partial steps, cf. Fig.\Qleft. 
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Figure 5: Evolution of spatial distribution of the state (u, 7T, £) with also the von Mises 
stress and the residuum R from (14.17ft at selected instants for the symmetric geometry from 
Fig.{2\- right. The deformation is visualized by a displacement u magnified 250 x, and t — 0.1 
was used. The process inherits the symmetry of the specimen and loading. In contrast to 
Fig.0 damage occurs rather later on the left-bottom corner and propagates fast, hence the 
snapshots are selected not in an equidistant way here. 
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Figure 6: Evolution of the dissipated energy Diss^(z; [0, t]) (top 3 curves) and integrated 
residua f* f n Rdxdt (bottom 3 curves) in two experiments from Figured In particular, it 
again shows good tendency of convergence and, moreover, that the violation of the approx¬ 
imate maximum-dissipation principle is small with respect to the overall dissipated energy 
and the evolution was stress driven with a good accuracy about 1-2%. 
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-time step=10~ 2 
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Figure 7: A detailed scaling of the bottom 3 curves (—the residua in AMDP) from Figure 
A good convergence to zero is seen in plastihcation period while some residuum is generated 
during damage period where nonconvexity of the stored energy truly comes into effect. 


it may be that there is no solution inheritting this symmetry or that experimental evidence 
shows preferences for nonsymmetric solutions. Cf. the discussion in [T2, i.26]. In addition, the 
geometry in Fig. (2J- left would lead to a 2 x 1 rectangle with a partial “cut” in the mid-bottom 
side, which is not a Lipschitz domain. 

Remark 6.2 (Recovery of the integrated maximum-dissipation principle IMDP). It should 
be emphasized that, even if the intuitively straightforward AMDP is asymptotically satisfied, 
the recovery of even the less-selective IMDP (13.811 for r —> 0 is not clear. This is obviously 
related with instability of IMDP under data perturbation if $(t, •) is not convex. Here, to 
recover the IMDP on I, it would suffice to show that for all £ > 0 there is r £ > 0 such that 
for any 0 < r < r e it holds 

T/r 

inf (£(t), z(kr)—z(kr— r)) — t%(z(kr)—z(kr— r)) > — e (6.1) 
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for some selection £(£) G —d z ^(t,u(t), z(t)), cf. the definition (13.71) and realize that the equi¬ 
distant partitions are cofinal in all partitions of I. This can be guaranteed only under rather 
strong conditions, namely if, for all £ > 0, there is r £ > 0 such that for any 0 < r < r e , the 
following strengthened version of the AMDP 

T/t £ 

YMMtk) - Z T (tk- 1)) - (fr (t),Z T (tk) - Z T (tk-1 )) < £ (6.2) 

k =1 


holds for tk = kr e , any tk-i < t < tk, and some £ T (t) G —d z £’(t,u T (t),z T (t)), and if 
£ T (t) —*■£(£), which can be assumed due to the available a-priori estimates used in the proof of 
Proposition 14.21 Using also (14.71) and the (norm,weak)-upper semicontinuity of d z S(t , •, •), in 
the limit for r —> 0, from such a strengthened AMDP, one can read Yk=i &( z {tk)—z(tk~ i)) — 
(£(£), z(tk) — z(tk- 1 )) < £ for any tk -i < t < tk and for some £(£) G —d z S(t,u[t), z(t)), from 
which (16. 1 j) indeed follows. In fact, our intuitive version of AMDP from Remark [O] computa¬ 
tionally verified (j6.2[) in Figure [7] in particular examples for t = r £ only. And, on top of it, we 
would need (16.2[) to be shown rather for G (—d 7r S > (t, u T (t), z T (t—r)), —d^Sit, u T (t), z T (t))). 
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